﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary.Statistics.RandomNumberGenerator;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the numerical integration of multivariate integrals by using 
    /// a Monte Carlo method.
    /// </summary>
    [Serializable]
    public class SimpleMultivariateMCIntegrator
    {
        /// <summary>
        /// This field holds the integral.
        /// </summary>
        private MultivariateIntegral integral;

        /// <summary>
        /// Initializes a new instance of the <see cref="SimpleMultivariateMCIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral to integrate.</param>
        public SimpleMultivariateMCIntegrator(MultivariateIntegral integral)
        {
            if (integral == (MultivariateIntegral) null)
            {
                throw new ArgumentNullException("integral");
            }

            if (integral.A == (GeneralVector) null)
            {
                throw new ArgumentNullException("integral");
            }

            if (integral.B == (GeneralVector) null)
            {
                throw new ArgumentNullException("integral");
            }

            this.integral = integral;
        }

        /// <summary>
        /// Gets or sets the integral to integrate.
        /// </summary>
        /// <value>The integral to integrate.</value>
        public MultivariateIntegral Integral
        {
            get { return this.integral; }
            set { this.integral = value; }
        }

        /// <summary>
        /// Integrates the specified integral in a numerical way.
        /// </summary>
        /// <param name="numberOfPoints">The maximum number of points to use for the integration.</param>
        /// <returns>The value of the integral.</returns>
        public double Integrate(int numberOfPoints)
        {
            if (this.integral.A.Count != this.integral.B.Count)
            {
                throw new ArgumentException("The number of integral border items a and b have to be even.");
            }

            double sum = 0.0;
            int dimension = this.integral.A.Count;
            GeneralVector point = new GeneralVector(dimension);
            XORShiftRandomGenerator random = new XORShiftRandomGenerator();

            for (int i = 0; i < numberOfPoints; i++)
            {
                for (int j = 0; j < dimension; j++)
                {
                    point[j] = this.integral.A[j] + random.NextDouble()*(this.integral.B[j] - this.integral.A[j]);
                }

                double fv = this.integral.Function.SolveAt(point);

                sum += fv;
            }

            double volume = 1.0;

            for (int i = 0; i < dimension; i++)
            {
                volume *= (this.integral.B[i] - this.integral.A[i]);
            }

            return (volume*sum)/numberOfPoints;
        }
    }
}